##################################################################################################################
# Functions that implement Individual testing, Master pool testing, Dorfman testing, Array testing


#############################################################################
# R function: This function simulates individual level testing and stores the
#             testing responses in accordance to the data structure 
#             required to fit the Bayesian GT regression model presented
#             in McMahan et al. via the R function Bayes.reg.GT
#
# Input: 
# Y.true = The true statuses of the individuals
# Se = The testing sensitivity used for both pools and individual testing
# Sp = The testing specificity used for both pools and individual testing
# 

Ind.decode<-function(Y.true,Se,Sp){
N<-length(Y.true)
Y<-rbinom(N,1,(Se*Y.true+(1-Sp)*(1-Y.true)))
return(Y)
}



#####################################################################
# R function: This function simulates Initial pool testing and stores the 
#             testing responses in accordance to the data structure 
#             required to fit the Bayesian GT regression model presented
#             in McMahan et al. via the R function Bayes.reg.GT
#
# Input: 
# Y.true = The true statuses of the individuals
# Se = A vector of testing sensitivities, where the first element is the
#      testing sensitivity for the pools and the second entry is the 
#      test sensitivity for individual testing
# Sp = A vector of testing specificities, where the first element is the
#      testing specificity for the pools and the second entry is the 
#      test specificity for individual testing
# cj = pool size to be used (Note: The number of individuals N should be 
#      evenly divisible by cj, this is only for decoding purposes; i.e., 
#      the regression methods do not require this condition)


Pool.test<-function(Y.true,Se,Sp,cj){
N<-length(Y.true)
Jmax<-N/cj 
J<-1

Y<-matrix(-99,nrow=N, ncol=3) 
Z<-matrix(-99,nrow=Jmax,ncol=cj+3) 


for(j in 1:(N/cj)){
prob<-ifelse(sum(Y.true[((j-1)*cj+1):(j*cj)])>0,Se,1-Sp)
Z[J,1]<-rbinom(n=1,size=1,prob=prob)
Z[J,2]<-cj
Z[J,3]<-1
Z[J,4:(cj+3)]<-((j-1)*cj+1):(j*cj)
Y[((j-1)*cj+1):(j*cj),2]<-1
Y[((j-1)*cj+1):(j*cj),3]<-J
J<-J+1
}

J<-J-1
Z<-Z[1:J,]

return(list("Z"=Z, "Y"=Y))
}



#####################################################################
# R function: This function simulates Dorfman decoding and stores the 
#             testing responses in accordance to the data structure 
#             required to fit the Bayesian GT regression model presented
#             in McMahan et al. via the R function Bayes.reg.GT
#
# Input: 
# Y.true = The true statuses of the individuals
# Se = A vector of testing sensitivities, where the first element is the
#      testing sensitivity for the pools and the second entry is the 
#      test sensitivity for individual testing
# Sp = A vector of testing specificities, where the first element is the
#      testing specificity for the pools and the second entry is the 
#      test specificity for individual testing
# cj = pool size to be used (Note: The number of individuals N should be 
#      evenly divisible by cj, this is only for decoding purposes; i.e., 
#      the regression methods do not require this condition)


Dorfman.decode.diff.error<-function(Y.true,Se,Sp,cj){
N<-length(Y.true)
Jmax<-N+N/cj 
J<-1

Y<-matrix(-99,nrow=N, ncol=4) 
Z<-matrix(-99,nrow=Jmax,ncol=cj+3) 


for(j in 1:(N/cj)){
prob<-ifelse(sum(Y.true[((j-1)*cj+1):(j*cj)])>0,Se[1],1-Sp[1])
Z[J,1]<-rbinom(n=1,size=1,prob=prob)
Z[J,2]<-cj
Z[J,3]<-1
Z[J,4:(cj+3)]<-((j-1)*cj+1):(j*cj)
Y[((j-1)*cj+1):(j*cj),2]<-1
Y[((j-1)*cj+1):(j*cj),3]<-J
J<-J+1
if(Z[J-1,1]==1){
for(k in ((j-1)*cj+1):(j*cj)){
prob<-ifelse(Y.true[k]>0,Se[2],1-Sp[2])
Z[J,1]<- rbinom(n=1,size=1,prob=prob)
Z[J,2]<-1
Z[J,3]<-2
Z[J,4]<-k
Y[k,2]<-2
Y[k,4]<-J
J<-J+1
}
}
}

J<-J-1
Z<-Z[1:J,]

return(list("Z"=Z, "Y"=Y))
}


#####################################################################
# R function: This function simulates Dorfman decoding and stores the 
#             testing responses in accordance to the data structure 
#             required to fit the Bayesian GT regression model presented
#             in McMahan et al. via the R function Bayes.reg.GT
#
# Input: 
# Y.true = The true statuses of the individuals
# Se = The testing sensitivity used for both pools and individual testing
# Sp = The testing specificity used for both pools and individual testing
# cj = pool size to be used (Note: The number of individuals N should be 
#      evenly divisible by cj, this is only for decoding purposes; i.e., 
#      the regression methods do not require this condition)



Dorfman.decode.same.error<-function(Y.true,Se,Sp,cj){
N<-length(Y.true)
Jmax<-N+N/cj 
J<-1

Y<-matrix(-99,nrow=N, ncol=4) 
Z<-matrix(-99,nrow=Jmax,ncol=cj+3) 


for(j in 1:(N/cj)){
prob<-ifelse(sum(Y.true[((j-1)*cj+1):(j*cj)])>0,Se,1-Sp)
Z[J,1]<-rbinom(n=1,size=1,prob=prob)
Z[J,2]<-cj
Z[J,3]<-1
Z[J,4:(cj+3)]<-((j-1)*cj+1):(j*cj)
Y[((j-1)*cj+1):(j*cj),2]<-1
Y[((j-1)*cj+1):(j*cj),3]<-J
J<-J+1
if(Z[J-1,1]==1){
for(k in ((j-1)*cj+1):(j*cj)){
prob<-ifelse(Y.true[k]>0,Se,1-Sp)
Z[J,1]<- rbinom(n=1,size=1,prob=prob)
Z[J,2]<-1
Z[J,3]<-1
Z[J,4]<-k
Y[k,2]<-2
Y[k,4]<-J
J<-J+1
}
}
}

J<-J-1
Z<-Z[1:J,]

return(list("Z"=Z, "Y"=Y))
}


#####################################################################
# R function: This function simulates Array decoding and stores the 
#             testing responses in accordance to the data structure 
#             required to fit the Bayesian GT regression model presented
#             in McMahan et al. via the R function Bayes.reg.GT
#
# Input: 
# Y.true = The true statuses of the individuals
# Se = A vector of testing sensitivities, where the first element is the
#      testing sensitivity for the pools and the second entry is the 
#      test sensitivity for individual testing
# Sp = A vector of testing specificities, where the first element is the
#      testing specificity for the pools and the second entry is the 
#      test specificity for individual testing
# cj = row and column pool sizes to be used (Note: The number of individuals 
#      N should be evenly divisible by cj^2, this is only for decoding 
#      purposes; i.e., the regression methods do not require this condition)


Array.decode.diff.error<-function(Y.true, Se, Sp, cj){

N<-length(Y.true)
Jmax<-2*N/cj +N
J<-1
AT<-N/(cj^2)

Y<-matrix(-99,nrow=N, ncol=5) 
Z<-matrix(-99,nrow=Jmax,ncol=cj+3) 

Y.A<-array(-99,c(cj,cj,AT))
ID.A<-array(-99,c(cj,cj,AT))
ind<-1
for(i in 1:AT){
for(j in 1:cj){
for(m in 1:cj){
Y.A[m,j,i]<-Y.true[ind]
ID.A[m,j,i]<-ind
ind<-ind+1
}}}



for(s in 1:AT){

array.yk<-Y.A[,,s]
array.id<-ID.A[,,s]

a<-rep(0,nrow(array.yk))
b<-rep(0,ncol(array.yk))

for(i in 1:cj){
   prob<- ifelse(sum(array.yk[i,])>0, Se[1], 1-Sp[1])
   g<- rbinom(n=1,size=1,prob=prob)
   a[i]<-g
   Z[J,1]<-g 
   Z[J,2]<-cj 
   Z[J,3]<-1
   Z[J,4:(cj+3)]<-array.id[i,]
   Y[array.id[i,],2]<-2
   Y[array.id[i,],3]<-J
   J<-J+1
}
for(j in 1:cj){
   prob<- ifelse(sum(array.yk[,j])>0, Se[1], 1-Sp[1])
   g<- rbinom(n=1,size=1,prob=prob)
   b[j]<-g
   Z[J,1]<-g 
   Z[J,2]<-cj 
   Z[J,3]<-1
   Z[J,4:(cj+3)]<-array.id[,j]
   Y[array.id[,j],4]<-J
   J<-J+1
}


if(sum(a)>0 & sum(b)>0){
array.yk1<-as.matrix(array.yk[(a==1),(b==1)])
array.id1<-as.matrix(array.id[(a==1),(b==1)])
for(i in 1:nrow(array.yk1)){
for(j in 1:ncol(array.yk1)){
   prob<- ifelse(array.yk1[i,j]>0, Se[2], 1-Sp[2])
   g<- rbinom(n=1,size=1,prob=prob)
   Z[J,1]<-g 
   Z[J,2]<-1 
   Z[J,3]<-2
   Z[J,4]<-array.id1[i,j]
   Y[array.id1[i,j],2]<-3
   Y[array.id1[i,j],5]<-J
   J<-J+1
}}}



if(sum(a)>0 & sum(b)==0){
array.yk1<-as.matrix(array.yk[(a==1),])
array.id1<-as.matrix(array.id[(a==1),])
for(i in 1:nrow(array.yk1)){
for(j in 1:ncol(array.yk1)){
   prob<- ifelse(array.yk1[i,j]>0, Se[2], 1-Sp[2])
   g<- rbinom(n=1,size=1,prob=prob)
   Z[J,1]<-g 
   Z[J,2]<-1 
   Z[J,3]<-2
   Z[J,4]<-array.id1[i,j]
   Y[array.id1[i,j],2]<-3
   Y[array.id1[i,j],5]<-J
   J<-J+1
}}}

if(sum(a)==0 & sum(b)>0){
array.yk1<-as.matrix(array.yk[,(b==1)])
array.id1<-as.matrix(array.id[,(b==1)])
for(i in 1:nrow(array.yk1)){
for(j in 1:ncol(array.yk1)){
   prob<- ifelse(array.yk1[i,j]>0, Se[2], 1-Sp[2])
   g<- rbinom(n=1,size=1,prob=prob)
   Z[J,1]<-g 
   Z[J,2]<-1 
   Z[J,3]<-2
   Z[J,4]<-array.id1[i,j]
   Y[array.id1[i,j],2]<-3
   Y[array.id1[i,j],5]<-J
   J<-J+1
}}}

} 

J<-J-1
Z<-Z[1:J,]

return(list("Z"=Z, "Y"=Y))
} #End function



#####################################################################
# R function: This function simulates Array decoding and stores the 
#             testing responses in accordance to the data structure 
#             required to fit the Bayesian GT regression model presented
#             in McMahan et al. via the R function Bayes.reg.GT
#
# Input: 
# Y.true = The true statuses of the individuals
# Se = The testing sensitivity used for both pools and individual testing
# Sp = The testing specificity used for both pools and individual testing
# cj = row and column pool sizes to be used (Note: The number of individuals 
#      N should be evenly divisible by cj^2, this is only for decoding 
#      purposes; i.e., the regression methods do not require this condition)


Array.decode.same.error<-function(Y.true, Se, Sp, cj){

N<-length(Y.true)
Jmax<-2*N/cj +N
J<-1
AT<-N/(cj^2)

Y<-matrix(-99,nrow=N, ncol=5) 
Z<-matrix(-99,nrow=Jmax,ncol=cj+3) 

Y.A<-array(-99,c(cj,cj,AT))
ID.A<-array(-99,c(cj,cj,AT))
ind<-1
for(i in 1:AT){
for(j in 1:cj){
for(m in 1:cj){
Y.A[m,j,i]<-Y.true[ind]
ID.A[m,j,i]<-ind
ind<-ind+1
}}}



for(s in 1:AT){

array.yk<-Y.A[,,s]
array.id<-ID.A[,,s]

a<-rep(0,nrow(array.yk))
b<-rep(0,ncol(array.yk))

for(i in 1:cj){
   prob<- ifelse(sum(array.yk[i,])>0, Se, 1-Sp)
   g<- rbinom(n=1,size=1,prob=prob)
   a[i]<-g
   Z[J,1]<-g 
   Z[J,2]<-cj 
   Z[J,3]<-1
   Z[J,4:(cj+3)]<-array.id[i,]
   Y[array.id[i,],2]<-2
   Y[array.id[i,],3]<-J
   J<-J+1
}
for(j in 1:cj){
   prob<- ifelse(sum(array.yk[,j])>0, Se, 1-Sp)
   g<- rbinom(n=1,size=1,prob=prob)
   b[j]<-g
   Z[J,1]<-g 
   Z[J,2]<-cj 
   Z[J,3]<-1
   Z[J,4:(cj+3)]<-array.id[,j]
   Y[array.id[,j],4]<-J
   J<-J+1
}


if(sum(a)>0 & sum(b)>0){
array.yk1<-as.matrix(array.yk[(a==1),(b==1)])
array.id1<-as.matrix(array.id[(a==1),(b==1)])
for(i in 1:nrow(array.yk1)){
for(j in 1:ncol(array.yk1)){
   prob<- ifelse(array.yk1[i,j]>0, Se, 1-Sp)
   g<- rbinom(n=1,size=1,prob=prob)
   Z[J,1]<-g 
   Z[J,2]<-1 
   Z[J,3]<-1
   Z[J,4]<-array.id1[i,j]
   Y[array.id1[i,j],2]<-3
   Y[array.id1[i,j],5]<-J
   J<-J+1
}}}



if(sum(a)>0 & sum(b)==0){
array.yk1<-as.matrix(array.yk[(a==1),])
array.id1<-as.matrix(array.id[(a==1),])
for(i in 1:nrow(array.yk1)){
for(j in 1:ncol(array.yk1)){
   prob<- ifelse(array.yk1[i,j]>0, Se, 1-Sp)
   g<- rbinom(n=1,size=1,prob=prob)
   Z[J,1]<-g 
   Z[J,2]<-1 
   Z[J,3]<-1
   Z[J,4]<-array.id1[i,j]
   Y[array.id1[i,j],2]<-3
   Y[array.id1[i,j],5]<-J
   J<-J+1
}}}

if(sum(a)==0 & sum(b)>0){
array.yk1<-as.matrix(array.yk[,(b==1)])
array.id1<-as.matrix(array.id[,(b==1)])
for(i in 1:nrow(array.yk1)){
for(j in 1:ncol(array.yk1)){
   prob<- ifelse(array.yk1[i,j]>0, Se, 1-Sp)
   g<- rbinom(n=1,size=1,prob=prob)
   Z[J,1]<-g 
   Z[J,2]<-1 
   Z[J,3]<-1
   Z[J,4]<-array.id1[i,j]
   Y[array.id1[i,j],2]<-3
   Y[array.id1[i,j],5]<-J
   J<-J+1
}}}

} 

J<-J-1
Z<-Z[1:J,]

return(list("Z"=Z, "Y"=Y))
} #End function


###########################################
# Function used to compute the coverage
# probabilities
cov.func<-function(data,sim,P,alpha,true){
cov<-matrix(-99,sim,P)
for(s in 1:sims){
for(p in 1:P){
int<-quantile(data[,p,s],c(alpha/2,1-alpha/2))
cov[s,p]<-int[1]<=true[p] & true[p]<=int[2]
}}
cov<-apply(cov,2,mean)
return(cov)
}